!======================================================
!                      v0(j),var(j)
! ------------#------------ j ------------#------------
!      fm(j-1)|fp(j-1)               fm(j)|fp(j)
!           flux(j-1)                 flux(j)
!
! fp(j): ---(-2)---(-1)---(0)---(1)---(2)---(3)---
!         1: |--------------|
!                 2: |-----------|
!                       3: |-----------|
!
! fm(j): ---(-2)---(-1)---(0)---(1)---(2)---(3)---
!                             1: |-----------|
!                       2: |-----------|
!                 3: |-----------|
!======================================================


subroutine wenozn5_p( fps, fpx )
    use commondata
    implicit none
    
    integer :: i
    real*8 :: fps(-2:2), fpx
    real*8 :: omega(3), fu(3), beta(3)
    real*8 :: tao, tao5, tao8, ac, aa
    real*8 :: totalpha, alpha(3), fomega(3)
    
    omega(1) = 0.1d0;   omega(2) = 0.6d0;   omega(3) = 0.3d0
    
    fu(1) =  1.0d0/3.0d0*fps(-2) - 7.0d0/6.0d0*fps(-1) + 11.0d0/6.0d0*fps(0)
    fu(2) = -1.0d0/6.0d0*fps(-1) + 5.0d0/6.0d0*fps(0)  + 1.0d0/3.0d0*fps(1)
    fu(3) =  1.0d0/3.0d0*fps(0)  + 5.0d0/6.0d0*fps(1)  - 1.0d0/6.0d0*fps(2)
    
    beta(1) = 13.0d0/12.0d0*( fps(-2) - 2.0d0*fps(-1) + fps(0) )**2 &
            + 0.25d0*( fps(-2) - 4.0d0*fps(-1) + 3.0d0*fps(0) )**2
    beta(2) = 13.0d0/12.0d0*( fps(-1) - 2.0d0*fps(0) + fps(1) )**2 &
            + 0.25d0*( fps(1) - fps(-1) )**2
    beta(3) = 13.0d0/12.0d0*( fps(0) - 2.0d0*fps(1) + fps(2) )**2 &
            + 0.25d0*( 3.0d0*fps(0) - 4.0d0*fps(1) + fps(2) )**2
    
    tao5 = abs(beta(1) - beta(3))
    tao = tao5
    aa = 10.d0
    ac = aa*( ( beta(1) + beta(3) - tao + eps )/( tao + eps ) )**2
    tao8 = ( fps(-2) - 4.d0*fps(-1) + 6.d0*fps(0) - 4.d0*fps(1) + fps(2) )**2
    
    alpha(:) = omega(:)*( ac + tao8/(beta(:)+eps) )
    totalpha = alpha(1) + alpha(2) + alpha(3)
    fomega(:) = alpha(:)/totalpha
    fpx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)
    
    return
    end subroutine wenozn5_p


subroutine wenozn5_m( fms, fmx )
    use commondata
    implicit none
    
    integer :: i
    real*8 :: fms(-1:3), fmx
    real*8 :: omega(3), fu(3), beta(3)
    real*8 :: tao, tao5, tao8, ac, aa
    real*8 :: totalpha, alpha(3), fomega(3)
    
    omega(1) = 0.1d0;   omega(2) = 0.6d0;   omega(3) = 0.3d0
    
    fu(1) =  1.0d0/3.0d0*fms(3) - 7.0d0/6.0d0*fms(2) + 11.0d0/6.0d0*fms(1)
    fu(2) = -1.0d0/6.0d0*fms(2) + 5.0d0/6.0d0*fms(1) + 1.0d0/3.0d0*fms(0)
    fu(3) =  1.0d0/3.0d0*fms(1) + 5.0d0/6.0d0*fms(0) - 1.0d0/6.0d0*fms(-1)
    
    beta(1) = 13.0d0/12.0d0*( fms(3) - 2.0d0*fms(2) + fms(1) )**2 &
            + 0.25d0*( fms(3) - 4.0d0*fms(2) + 3.0d0*fms(1) )**2
    beta(2) = 13.0d0/12.0d0*( fms(2) - 2.0d0*fms(1) + fms(0) )**2 &
            + 0.25d0*( fms(2) - fms(0) )**2
    beta(3) = 13.0d0/12.0d0*( fms(1) - 2.0d0*fms(0) + fms(-1) )**2 &
            + 0.25d0*( 3.0d0*fms(1) - 4.0d0*fms(0) + fms(-1) )**2
    
    tao5 = abs(beta(1) - beta(3))
    tao = tao5
    aa = 10.d0
    ac = aa*( ( beta(1) + beta(3) - tao + eps )/( tao + eps ) )**2
    tao8 = ( fms(-1) - 4.d0*fms(0) + 6.d0*fms(1) - 4.d0*fms(2) + fms(3) )**2
    
    alpha(:) = omega(:)*( ac + tao8/(beta(:)+eps) )
    totalpha = alpha(1) + alpha(2) + alpha(3)
    fomega(:) = alpha(:)/totalpha
    fmx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)
    
    return
    end subroutine wenozn5_m